%% produces the transition dynamics of the economy between steady states

A_s_new    = A_vec_s(ip,ip2);
A_n_new    = A_vec_n(ip,ip2);

%% find the avg. productivity of surviving firms at entry
f_entry_s = zeros(Ntot,1);
f_entry_n = zeros(Ntot,1);
f_entry_s(inde_bi_s)  = (speye(length(inde_bi_s))+ A_T_s(inde_bi_s,inde_bi_s)*dt)\h_00_s(inde_bi_s);%mass of mature firms per worker
f_entry_n(inde_bi_n)  = (speye(length(inde_bi_n))+ A_T_n(inde_bi_n,inde_bi_n)*dt)\h_00_n(inde_bi_n);%mass of mature firms per worker
Ez_entry_s            = sum(Z_s(inde_bi_s).*f_entry_s(inde_bi_s));
Ez_entry_n            = sum(Z_n(inde_bi_n).*f_entry_n(inde_bi_n));

%% compute output and average firm productivity at t=0
if demand_ext==1
    Y_new      =  (wei_s*l_s*A_s_new^(1/(1-nu)) + wei_n*l_n*A_n_new^(1/(1-nu)))^((nu-1)/(nu-2));
else
    Y_new      = 1;
end
w_s        = (nu-1)/nu*(A_s_new/Y_new)^(1/(1-nu));
w_n        = (nu-1)/nu*(A_n_new/Y_new)^(1/(1-nu));
Ez_s       = (w_s/(1-1/nu))^nu/Y_new;%avg quality per worker across all firms
Ez_n       = (w_n/(1-1/nu))^nu/Y_new;%avg quality per worker across all firms

%% compute the skake out and the mass of entrants on impact
f_m1_s     = f_dist_ss_s;
f_m1_n     = f_dist_ss_n;
f_dist_ss_s= f_m1_s;
f_dist_ss_n= f_m1_n;
get_shakout_distribution;
m_entry0_s = (Ez_s - sum(Z_s.*f_dist0_s))/(Ez_entry_s);
m_entry0_n = (Ez_n - sum(Z_n.*f_dist0_n))/(Ez_entry_n);
%here we are computing the shift in the distribution after the shock
if ip2==ip2_sub && ip ==ip1_sub  && incid_expe == 1
    f_plus_s     = (f_dist0_s + m_entry0_s*h_00_s);
    f_plus_n     = (f_dist0_n + m_entry0_n*h_00_n);
    f_0_s        = f_dist0_s;
    f_0_n        = f_dist0_n;
    get_figure_distribution;
end
%% simulations
Y_t    = ones(T_sim,1);
f_t_s  = zeros(Ntot,T_sim);
f_t_n  = zeros(Ntot,T_sim);

%notice: start initial distribution with SO effect
f_m1_s     = f_dist0_s*l_s;%mass of firms per province
f_m1_s(inde_bb_s)=0;
f_m1_n     = f_dist0_n*l_n;%mass of firms per province
f_m1_n(inde_bb_n)=0;

time=0;inde_sim=0;l_s=l_s_ss;l_n=l_n_ss;
while inde_sim<T_sim
    time=time+dt;inde_sim=inde_sim+1;

    %%%%%%%%%%%%%%%%%%%%%%% predetermined output
    y_s        = A_s_new*Ez_s;
    y_n        = A_n_new*Ez_n;
    chi_firm_s= chi_s*Ez_s;
    chi_firm_n= chi_n*Ez_n;

    %%%%%%%%%%%%%%%%%%%%%% distribution of non-bankruptcy/exogenous-exiting firms
    f_mature_s            = zeros(Ntot,1);
    f_mature_n            = zeros(Ntot,1);
    f_mature_s(inde_bi_s) = (speye(length(inde_bi_s))+ A_T_s(inde_bi_s,inde_bi_s)*dt)\f_m1_s(inde_bi_s);%mass of mature firms per worker
    f_mature_n(inde_bi_n) = (speye(length(inde_bi_n))+ A_T_n(inde_bi_n,inde_bi_n)*dt)\f_m1_n(inde_bi_n);%mass of mature firms per worker
    Ez_mature_s           = sum(Z_s.*f_mature_s)/l_s;%avg quality per worker across mature firms per worker
    Ez_mature_n           = sum(Z_n.*f_mature_n)/l_n;%avg quality per worker across mature firms per worker

    %compute entry and investment: notice that in the first period
    %it accounts for both the SO effect and elasping of one %period;
    m_entry_s  = (Ez_s - (Ez_mature_s))/(Ez_entry_s);%entry per worker
    m_entry_n  = (Ez_n - (Ez_mature_n))/(Ez_entry_n);%entry per worker
    i_s        = k_s_ben*m_entry_s;%investment per worker
    i_n        = k_n_ben*m_entry_n;%investment per worker

    %compute exit
    m_firms_s= sum(f_m1_s); %mass of firms from previous period
    m_firms_n= sum(f_m1_n); %mass of firms from previous period
    m_exit_s = sum(f_m1_s)-sum(f_mature_s)+l_s*m_entry_s*(1-sum(f_entry_s));  %mass of  exiting firms from previous period
    m_exit_n = sum(f_m1_n)-sum(f_mature_n)+l_n*m_entry_n*(1-sum(f_entry_n));  %mass of  exiting firms from previous period

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% update distribution
    f_t_s(inde_bi_s,inde_sim)= f_mature_s(inde_bi_s) + f_entry_s(inde_bi_s)*m_entry_s*l_s;%mass of firms
    f_t_n(inde_bi_n,inde_sim)= f_mature_n(inde_bi_n) + f_entry_n(inde_bi_n)*m_entry_n*l_n;%mass of firms
    %investment per worker
    inve_s     = A_s_new*sum(Z_s.*inv_s.*f_t_s(:,inde_sim))/l_s;
    inve_n     = A_n_new*sum(Z_n.*inv_n.*f_t_n(:,inde_sim))/l_n;

    f_m1_s=zeros(Ntot,1);
    f_m1_n=zeros(Ntot,1);
    f_m1_s(inde_bi_s) = f_t_s(inde_bi_s,inde_sim);
    f_m1_n(inde_bi_n) = f_t_n(inde_bi_n,inde_sim);

    %Save results

    welf_t_all(inde_sim,ip,ip2)  = (y_s-i_s-chi_firm_s-h_s*l_s^(eta)*dt-inve_s)*l_s*wei_s+(y_n-i_n-chi_firm_n-h_n*l_n^(eta)*dt-inve_n)*l_n*wei_n;

    if  inde_dilution==1 && l_reallocation==1 &&  demand_ext==1 && shiftdist==1 
        sub_time_s_mat(inde_sim,ip,ip2)   = i_s*size_subsidy*(1-sub_tax_s);
        sub_time_n_mat(inde_sim,ip,ip2)   = i_s*size_subsidy*(1-sub_tax_s)*(l_s*wei_s)/(l_n*wei_n);
        sub_time_mat(inde_sim,ip,ip2)     = sub_time_s_mat(inde_sim,ip,ip2)*l_s*wei_s-sub_time_n_mat(inde_sim,ip,ip2)*l_n*wei_n;
        W_time_s_vec(inde_sim,ip,ip2)     = y_s-i_s-chi_firm_s-h_s*l_s^(eta)*dt-inve_s+sub_time_s_mat(inde_sim,ip,ip2);
        W_time_n_vec(inde_sim,ip,ip2)     = y_n-i_n-chi_firm_n-h_n*l_n^(eta)*dt-inve_n-sub_time_n_mat(inde_sim,ip,ip2);
        W_time_vec(inde_sim,ip,ip2)       = W_time_s_vec(inde_sim,ip,ip2)*l_s*wei_s+W_time_n_vec(inde_sim,ip,ip2)*l_n*wei_n;
        y_time_mat(inde_sim,ip,ip2)       = y_s*l_s*wei_s+y_n*l_n*wei_n;
        y_time_s_mat(inde_sim,ip,ip2)     = y_s;
        y_time_n_mat(inde_sim,ip,ip2)     = y_n;
        i_time_mat(inde_sim,ip,ip2)       = (inve_s+i_s)*l_s*wei_s+(inve_n+i_n)*l_n*wei_n;
        i_time_s_mat(inde_sim,ip,ip2)     = inve_s+i_s;
        i_time_n_mat(inde_sim,ip,ip2)     = inve_n+i_n;
        inve_time_s_mat(inde_sim,ip,ip2)  = inve_s;
        inve_time_n_mat(inde_sim,ip,ip2)  = inve_n;
        chi_time_mat(inde_sim,ip,ip2)    = chi_firm_s*l_s*wei_s+chi_firm_n*l_n*wei_n;
        chi_time_s_mat(inde_sim,ip,ip2)  = chi_firm_s;
        chi_time_n_mat(inde_sim,ip,ip2)  = chi_firm_n;
        conj_time_mat(inde_sim,ip,ip2)    = h_n*l_n^(eta)*dt*l_n*wei_n+h_s*l_s^(eta)*dt*l_s*wei_s;
        conj_time_n_mat(inde_sim,ip,ip2)  = h_n*l_n^(eta)*dt;
        conj_time_s_mat(inde_sim,ip,ip2)  = h_s*l_s^(eta)*dt;
        l_time_s_mat(inde_sim,ip,ip2)     = l_s;
        l_time_n_mat(inde_sim,ip,ip2)     = l_n;
        entry_t_s(inde_sim,ip,ip2)        = m_entry_s*l_s/m_firms_s;  %entry rate per firm
        entry_t_n(inde_sim,ip,ip2)        = m_entry_n*l_n/m_firms_n;  %entry rate per firm
        exit_t_s(inde_sim,ip,ip2)         = m_exit_s/sum(f_m1_s); %exit rate per firm
        exit_t_n(inde_sim,ip,ip2)         = m_exit_n/sum(f_m1_n); %exit rate per firm
        debt0_t_s(inde_sim,ip,ip2)        = debt0_va_s_vec(ip,ip2);
        debt0_t_n(inde_sim,ip,ip2)        = debt0_va_n_vec(ip,ip2);
        debt_t_n(inde_sim,ip,ip2)         = bb_vec_mat_n'*f_m1_n/sum(f_m1_n)*A_n_new/R_vec_n(ip,ip2);
        debt_t_s(inde_sim,ip,ip2)         = bb_vec_mat_s'*f_m1_s/sum(f_m1_s)*A_s_new/R_vec_s(ip,ip2);
    end

    if ip2==1 &&  ip==1 && inde_sim==1 
        cons_ss=(y_s-i_s-inve_s)*l_s*wei_s+(y_n-i_n-inve_n)*l_n*wei_n;
    end

    %update labor
    if l_reallocation==1
        dl_s     = l_s_vec(ip,ip2)-l_s_ss;
        l_s_new  = l_s_ss + dl_s*(1-exp(-psi*time));
        l_n_new  = (1-l_s_new*wei_s)/wei_n;
        l_n      = l_n_new;
        l_s      = l_s_new;
    end

    %update output and average productovity
    if demand_ext==1
        Y_new      = (wei_s*l_s*A_s_new^(1/(1-nu)) + wei_n*l_n*A_n_new^(1/(1-nu)))^((nu-1)/(nu-2));
    else
        Y_new      = 1;
    end
    w_s        = (nu-1)/nu*(A_s_new/Y_new)^(1/(1-nu));
    w_n        = (nu-1)/nu*(A_n_new/Y_new)^(1/(1-nu));
    Ez_s       = (w_s/(1-1/nu))^nu/Y_new;%avg quality per worker across all firms
    Ez_n       = (w_n/(1-1/nu))^nu/Y_new;%avg quality per worker across all firms
end
debt0_s(ip,ip2)        = debt0_va_s_vec(ip,ip2);
debt0_n(ip,ip2)        = debt0_va_n_vec(ip,ip2);
debt0t_s(ip,ip2)       = debt0t_va_s_vec(ip,ip2);

%% save results per parameter combination
disc                   = (1/(1+r)).^(0:dt:T_y-dt);
welf_disc_all(ip,ip2)  = r*((disc*welf_t_all(:,ip,ip2))+disc(end)/r*welf_t_all(end,ip,ip2));

% Save PDV of variables in baseline case
if  inde_dilution  ==1 && l_reallocation==1 &&  demand_ext==1 && shiftdist==1  && incid_expe == 1
    welf_disc(ip,ip2)      =r*((disc*W_time_vec(:,ip,ip2))+disc(end)/r*W_time_vec(end,ip,ip2));
    welf_disc_s(ip,ip2)    =r*((disc*W_time_s_vec(:,ip,ip2))+disc(end)/r*W_time_s_vec(end,ip,ip2));
    welf_disc_n(ip,ip2)    =r*((disc*W_time_n_vec(:,ip,ip2))+disc(end)/r*W_time_n_vec(end,ip,ip2));
    welf_disc_ss(ip,ip2)   =W_time_vec(end,ip,ip2);
    welf_disc_s_ss(ip,ip2) =W_time_s_vec(end,ip,ip2);
    welf_disc_n_ss(ip,ip2) =W_time_n_vec(end,ip,ip2);
    y_disc(ip,ip2)         =r*((disc*y_time_mat(:,ip,ip2))+disc(end)/r*y_time_mat(end,ip,ip2));
    y_disc_s(ip,ip2)       =r*((disc*y_time_s_mat(:,ip,ip2))+disc(end)/r*y_time_s_mat(end,ip,ip2));
    y_disc_n(ip,ip2)       =r*((disc*y_time_n_mat(:,ip,ip2))+disc(end)/r*y_time_n_mat(end,ip,ip2));
    i_disc(ip,ip2)         =r*((disc*i_time_mat(:,ip,ip2))+disc(end)/r*i_time_mat(end,ip,ip2));
    i_disc_s(ip,ip2)       =r*((disc*i_time_s_mat(:,ip,ip2))+disc(end)/r*i_time_s_mat(end,ip,ip2));
    i_disc_n(ip,ip2)       =r*((disc*i_time_n_mat(:,ip,ip2))+disc(end)/r*i_time_n_mat(end,ip,ip2));
    chi_disc(ip,ip2)      =r*((disc*chi_time_mat(:,ip,ip2))+disc(end)/r*chi_time_mat(end,ip,ip2));
    chi_disc_s(ip,ip2)    =r*((disc*chi_time_s_mat(:,ip,ip2))+disc(end)/r*chi_time_s_mat(end,ip,ip2));
    chi_disc_n(ip,ip2)    =r*((disc*chi_time_n_mat(:,ip,ip2))+disc(end)/r*chi_time_n_mat(end,ip,ip2));
    conj_disc(ip,ip2)      =r*((disc*conj_time_mat(:,ip,ip2))+disc(end)/r*conj_time_mat(end,ip,ip2));
    conj_disc_s(ip,ip2)    =r*((disc*conj_time_s_mat(:,ip,ip2))+disc(end)/r*conj_time_s_mat(end,ip,ip2));
    conj_disc_n(ip,ip2)    =r*((disc*conj_time_n_mat(:,ip,ip2))+disc(end)/r*conj_time_n_mat(end,ip,ip2));
    levratio_0(ip,ip2)     =debt0_va_s_vec(ip,ip2);
    sub_disc_s(ip,ip2)     =r*((disc*sub_time_s_mat(:,ip,ip2)+disc(end)/r*sub_time_s_mat(end,ip,ip2)));
    sub_disc_n(ip,ip2)     =r*((disc*sub_time_n_mat(:,ip,ip2)+disc(end)/r*sub_time_n_mat(end,ip,ip2)));
    sub_disc(ip,ip2)       =r*((disc*sub_time_mat(:,ip,ip2)+disc(end)/r*sub_time_mat(end,ip,ip2)));
    b0_s(ip,ip2)           = (zlow_s*b0_hig_s*q_hig_s+zhig_s*b0_low_s*q_low_s+zmid_s*b0_mid_s*q_mid_s)/(1/nu*A_s_new-chi_s);
end

%if abs((ie-4).*(ie-6))>0.01, ip=length(tau_vec);end
